//------------------------------------------------------------------
//  POLYEC.c
//------------------------------------------------------------------
/************************************************************************************
*																					*
*		Polynomial version of elliptic curve functions.  essentially the same as	*
*  optimal normal basis but uses polynomial functions.  Biggest difference is that	*
*  once poly_prime is picked (see polymain.c) you need to compute the "T matrix" to	*
*  help embed data onto a curve by solving a quadratic equation.					*
*																					*
************************************************************************************/
//-------------------------------------------------------------------
#include <stdio.h>
//-------------------------------------------------------------------
#include "Top.h"
#include "Bigint.h"
#include "BigintEC.h"
#include "Poly.h"
#include "PolyEC.h"
//-------------------------------------------------------------------
extern	FIELD2N poly_prime;
//------------------------------------------------------------------
//  Global parameters for polynomial math.  poly_prime is the irreducible polynomial
//	for the Galois field arithmetic.  Tmatrix is used to solve quadratic equations
//	and Smatrix is used to compute square roots.
//------------------------------------------------------------------
/****************************************************************************
*   Implement elliptic curve point addition for polynomial basis form.  	*
*  This follows R. Schroeppel, H. Orman, S. O'Mally, "Fast Key Exchange with*
*  Elliptic Curve Systems", CRYPTO '95, TR-95-03, Univ. of Arizona, Comp.   *
*  Science Dept.                                                            *
****************************************************************************/
void poly_eneg(POLY_POINT *p)
{
	register INDEX i;
	
	SUMLOOP(i) p->y.e[i] ^= p->x.e[i];
}
//------------------------------------------------------------------
//
//------------------------------------------------------------------
void poly_esum(POLY_POINT* p1, POLY_POINT* p2, POLY_POINT* p3, FIELD2N* irr_poly, POLY_CURVE* curv)
{
    INDEX   i;
    FIELD2N  x1, y1, theta, onex, theta2;
    ELEMENT  check;
	
/*  check if p1 or p2 is point at infinity  */

	check = 0;
	SUMLOOP(i) check |= p1->x.e[i] | p1->y.e[i];
	if (!check)
	{
		printf("summing point at infinity\n");
		poly_copy_point( p2, p3);
		return;
	}
	check = 0;
	SUMLOOP(i) check |= p2->x.e[i] | p2->y.e[i];
	if (!check) 
	{
		printf("summing point at infinity\n");
		poly_copy_point( p1, p3);
		return;
	}
	
/*  compute theta = (y_1 + y_2)/(x_1 + x_2)  */

    poly_null(&x1);
    poly_null(&y1);
    check = 0;
    SUMLOOP(i) 
    {
		x1.e[i] = p1->x.e[i] ^ p2->x.e[i];
		y1.e[i] = p1->y.e[i] ^ p2->y.e[i];
		check |= x1.e[i];
    }
    if (!check)   /* return point at infinity  */
    {
    	printf("input to elliptic sum has duplicate x values\n"); 
    	poly_null(&p3->x);
    	poly_null(&p3->y);
    	return;
    }
    poly_mod_inv( &x1, &onex, irr_poly);
    poly_mod_mul( &onex, &y1, &theta, irr_poly);  /*  compute y_1/x_1 = theta */
    poly_mod_mul(&theta, &theta, &theta2, irr_poly);   /* then theta^2  */

/*  with theta and theta^2, compute x_3  */

    if (curv->form)
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i)
	    	p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ x1.e[i];

/*  next find y_3  */

    SUMLOOP (i) x1.e[i] = p1->x.e[i] ^ p3->x.e[i];
    poly_mod_mul( &x1, &theta, &theta2, irr_poly);
    SUMLOOP (i) p3->y.e[i] = theta2.e[i] ^ p3->x.e[i] ^ p1->y.e[i];
}
//------------------------------------------------------------------
//  elliptic curve doubling routine for Schroeppel's algorithm over polymomial
//  basis.  Enter with p1, p3 as source and destination as well as curv
//  to operate on.  Returns p3 = 2*p1.
//------------------------------------------------------------------
void poly_edbl(POLY_POINT* p1, POLY_POINT* p3, FIELD2N* irr_poly, POLY_CURVE* curv)
{
    FIELD2N  x1, y1, theta, theta2, t1;
    INDEX   i;
    ELEMENT  check;
    
    check = 0;
    SUMLOOP (i) check |= p1->x.e[i];
    if (!check)
    {
    	printf("doubling point at infinity\n");
    	poly_null(&p3->x);
    	poly_null(&p3->y);
    	return;
    }

/*  first compute theta = x + y/x  */

    poly_mod_inv( &p1->x, &x1, irr_poly);
    poly_mod_mul( &x1, &p1->y, &y1, irr_poly);
    SUMLOOP (i) theta.e[i] = p1->x.e[i] ^ y1.e[i];

/*  next compute x_3  */

    poly_mod_mul( &theta, &theta, &theta2, irr_poly);
    if(curv->form)
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i] ^ curv->a2.e[i];
    else
		SUMLOOP (i) p3->x.e[i] = theta.e[i] ^ theta2.e[i];

/*  and lastly y_3  */

    theta.e[NUMWORD] ^= 1;			/*  theta + 1 */
    poly_mod_mul( &theta, &p3->x, &t1, irr_poly);
    poly_mod_mul( &p1->x, &p1->x, &x1, irr_poly);
    SUMLOOP (i) p3->y.e[i] = x1.e[i] ^ t1.e[i];
}
//------------------------------------------------------------------
//  subtract two points on a curve.  just negates p2 and does a sum.
//  Returns p3 = p1 - p2 over curv.
//------------------------------------------------------------------
void poly_esub(POLY_POINT* p1, POLY_POINT* p2, POLY_POINT* p3, FIELD2N* irr_poly, POLY_CURVE* curv)
{
    POLY_POINT   negp;
    INDEX   i;

    poly_copy ( &p2->x, &negp.x);
    poly_null (&negp.y);
    SUMLOOP(i) negp.y.e[i] = p2->x.e[i] ^ p2->y.e[i];
    poly_esum (p1, &negp, p3, irr_poly, curv);
}
//------------------------------------------------------------------
//  need to move points around, not just values.  Optimize later.  
//------------------------------------------------------------------
void poly_copy_point(POLY_POINT* p1, POLY_POINT* p2)
{
	poly_copy (&p1->x, &p2->x);
	poly_copy (&p1->y, &p2->y);
}
//------------------------------------------------------------------
//  Routine to compute kP where k is an integer (base 2, not normal basis)
//	and P is a point on an elliptic curve.  This routine assumes that K
//	is representable in the same bit field as x, y or z values of P.  Since
//	the field size determines the largest possible order, this makes sense.
//  Enter with: integer k, source point P, curve to compute over (curv) 
//  Returns with: result point R.
//
//  Reference: Koblitz, "CM-Curves with good Cryptografic Properties", 
//	Springer-Verlag LNCS #576, p279 (pg 284 really), 1992
//------------------------------------------------------------------
void  poly_emul(FIELD2N* k, POLY_POINT* p, POLY_POINT* r, FIELD2N* irr_poly, POLY_CURVE* curv)
{
	char		blncd[NUMBITS+1];
	INDEX		bit_count, i;
	ELEMENT		notzero;
	FIELD2N		number;
	POLY_POINT		temp;

/*  make sure input multiplier k is not zero.
	Return point at infinity if it is.
*/
	poly_copy( k, &number);
	notzero = 0;
	SUMLOOP (i) notzero |= number.e[i];
	if (!notzero)
	{
		poly_null (&r->x);
		poly_null (&r->y);
		return;
	}

/*  convert integer k (number) to balanced representation.
	Called non-adjacent form in "An Improved Algorithm for
	Arithmetic on a Family of Elliptic Curves", J. Solinas
	CRYPTO '97. This follows algorithm 2 in that paper.
*/
	bit_count = 0;
	while (notzero)
	{
	/*  if number odd, create 1 or -1 from last 2 bits  */
	
		if ( number.e[NUMWORD] & 1 )
		{
			blncd[bit_count] = 2 - (number.e[NUMWORD] & 3);
			
	/*  if -1, then add 1 and propagate carry if needed  */
			
			if ( blncd[bit_count] < 0 )
			{
				for (i=NUMWORD; i>=0; i--)
				{
					number.e[i]++;
					if (number.e[i]) break;
				}
			}
		}
		else
			blncd[bit_count] = 0;
	
	/*  divide number by 2, increment bit counter, and see if done  */
	
		number.e[NUMWORD] &= ~0 << 1;
		poly_rot_right( &number);
		bit_count++;
		notzero = 0;
		SUMLOOP (i) notzero |= number.e[i];
	}
		
/*  now follow balanced representation and compute kP  */

	bit_count--;
	poly_copy_point(p,r);		/* first bit always set */
	while (bit_count > 0) 
	{
	  poly_edbl(r, &temp, irr_poly, curv);
	  bit_count--;
	  switch (blncd[bit_count]) 
	  {
	     case 1: poly_esum (p, &temp, r, irr_poly, curv);
				 break;
	     case -1: poly_esub (&temp, p, r, irr_poly, curv);
				  break;
	     case 0: poly_copy_point (&temp, r);
	   }
	}
}
//------------------------------------------------------------------
//
//------------------------------------------------------------------
void print_poly_field(char* string, FIELD2N* x)
{
	INDEX i;
	
	printf("%s\n",string);
	SUMLOOP(i) printf("%8x ",x->e[i]);
	printf("\n");
}
//------------------------------------------------------------------
//
//------------------------------------------------------------------
void print_poly_point(char* string, POLY_POINT* point)
{
	INDEX i;
	
	printf("%s\n",string);
	printf("x: ");
	SUMLOOP(i) printf("%8x ",point->x.e[i]);
	printf("\n");
	printf("y: ");
	SUMLOOP(i) printf("%8x ",point->y.e[i]);
	printf("\n");
}
//------------------------------------------------------------------
//
//------------------------------------------------------------------
void print_poly_curve(char* string, POLY_CURVE* curv)
{
	INDEX i;
	
	printf("%s\n",string);
	printf("form: %d\n",curv->form);
	if (curv->form)
	{
		printf("a2: ");
		SUMLOOP(i) printf("%lx ",curv->a2.e[i]);
		printf("\n");
	}
	printf("a6: ");
	SUMLOOP(i) printf("%lx ",curv->a6.e[i]);
	printf("\n\n");
}
//------------------------------------------------------------------
//
//------------------------------------------------------------------
void poly_esum_LD(FIELD2N* bpx, POLY_POINT_LD* P1, POLY_POINT_LD* P2, POLY_POINT_LD* P3)
{
	POLY_POINT_LD  psum;
	FIELD2N        temp1;
	FIELD2N        temp2;
	FIELD2N        temp3;
	FIELD2N        temp4;
	// Z3 = (X1*Z2 + X2*Z1)^2
	poly_mod_mul(&P1->X, &P2->Z,  &temp1,  &poly_prime); // temp1  = X1*Z2
	poly_mod_mul(&P2->X, &P1->Z,  &temp2,  &poly_prime); // temp2  = X2*Z1
	poly_mod_add(&temp1, &temp2,  &temp3,  &poly_prime); // temp3  = X1*Z2 + X2*Z1
	poly_mod_mul(&temp3, &temp3,  &psum.Z, &poly_prime); // psum.Z = (X1*Z2 + X2*Z1)^2
	// X3 = x*Z3 + (X1*Z2)*(X2*Z1)
	poly_mod_mul(&temp1, &temp2,  &temp3,  &poly_prime); // temp3  = (X1*Z2)*(X2*Z1)
	poly_mod_mul(bpx,    &psum.Z, &temp4,  &poly_prime); // temp4  = x*Z3
	poly_mod_add(&temp4, &temp3,  &psum.X, &poly_prime); // PSUM.x = x*Z3 + (X1*Z2)*(X2*Z1)
	//
	poly_copy(&psum.X, &P3->X);
	poly_copy(&psum.Z, &P3->Z);
}

void poly_edlb_LD(POLY_POINT_LD* P1, POLY_POINT_LD* P3, POLY_CURVE* curv)
{
	POLY_POINT_LD pdbl;
	FIELD2N       temp1;
	FIELD2N       temp2;
	FIELD2N       temp3;
	FIELD2N       temp4;
	FIELD2N       temp5;
	// Z3 = (X1*Z1)^2
	poly_mod_mul(&P1->X,    &P1->Z, &temp1,  &poly_prime);   // temp1  = X1*Z1
	poly_mod_mul(&temp1,    &temp1, &pdbl.Z, &poly_prime);   // pdbl.Z = (X1*Z1)^2
	// X3 = X1^4 + b*Z1^4
	poly_mod_mul(&P1->X,    &P1->X, &temp1,  &poly_prime);    // temp1  = X1^2
	poly_mod_mul(&temp1,    &temp1, &temp2,  &poly_prime);    // temp2  = X1^4
	poly_mod_mul(&P1->Z,    &P1->Z, &temp3,  &poly_prime);    // temp3  = Z1^2
	poly_mod_mul(&temp3,    &temp3, &temp4,  &poly_prime);    // temp4  = Z1^4
	poly_mod_mul(&curv->a6, &temp4, &temp5,  &poly_prime); // temp5  = b*Z1^4
	poly_mod_add(&temp2,    &temp5, &pdbl.X, &poly_prime);   // pdlb.X = X1^4 + b*Z1^4
	//
	poly_copy(&pdbl.X, &P3->X);
	poly_copy(&pdbl.Z, &P3->Z);
}